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Abstract 

We analytically construct an evolutionary potential landscape in the Wright- 
Fisher model. Compared to the classical fitness landscape, it gives a complete 
description for system's middle- and long- term behaviors. Compared to previous 
work, our formulation is valid in the whole parameter space and does not require 
a normalizable equilibrium distribution. We discuss the uphill and downhill dy- 
namics on the potential landscape and the associated timescales. When studying 
the average time to escape from an infinite potential peak, we find a new way to 
analytically approximate the result. Our results extend the use of Kramers' escape 
formula to the non-Gaussian cases and bridge previous results in two limits. We 
conclude that a divergent equilibrium distribution does not necessarily imply the 
fixation of an allele type. We show how further biological insights can be gained 
from our framework by discussing three related issues. 



1 Introduction 



The concept of adaptive landscape was first proposed by Wright (1932), in attempting to 



visualize his shifting balance theory in evolutionary dynamics. Since then, the metaphoric 
and visualizing part of this concept has been widely used in population genetics and 



evolutionary biology (Lande, 1976 Wright, 1988 Arnold et ah, 2001). Wright's original 



landscape can be interpreted as the "the surface of selective fitness" (Wright, 1967). 



Its definition has two main derivatives (Provine, 1986): project the individual fitness 



over space of allele/geno-/pheno- types or project the mean fitness of a population over 
the frequencies of allele/geno-/pheno- types. In this paper, we discuss the second type, 
and the landscape is defined over the allele frequency space. We do not project fitness, 
though, but define a measure of evolutionary potential; it includes the effects of other 
factors (mutation, migration, genetic drift, etc) besides selection. It does not coincide the 
classical definition of fitness landscape, but manifests the essential idea that is embedded 
in Wright's original proposition: There exist general laws underlying the evolutionary 
processes. 
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In terms of the middle- and long- term behaviors of an evolutionary system, fitness 
presents only part of the whole dynamical picture. Because of the non-trivial effect of 
random genetic drift, the classical fitness landscape does not necessarily coincide the equi- 
librium probability distribution with respect to their peak/valley numbers and locations. 
The mismatch exists even in the 1-D case and even if there is no other "directed" fac- 
tors (e.g. mutation/migration) besides selection. And when there is mutation/migration, 
finding a complete description for the evolutionary dynamics on the fitness landscape 



would be difficult (Gavrilets, 1997 Fear and Price, 1998) and sometimes artificial (Coyne 



et al. , 1997). An example is modeling the rare events and calculating its occurrence 



probability under the joint effects of mutation, selection and random drift (Lande, 1976 



Barton and Rouhani 1987; Orr 2009). In other aspects, lack of theoretical generality and 



self-consistency has put this classical concept under controversies (Pro vine, 1986 Rice 



2004 Kaplan, 2008). Even Wright (1988) himself claimed that people "was looking for 



something more mathematical than was intended (from the surfaces of selective value)". 

Despite above critiques on the fitness landscape, the heuristic usefulness and possible 
biological insights of certain scalar functions in evolution have gained more and more 
attentions. There are re-emerging interests in studying the evolutionary systems by 



constructing quantitative landscapes in specific contexts (Schluter and Nychka 
Waxman and Gavriletsj [2005| [Qianl [2005j [Hendry et al.j [2006110 arneiro and Hartll 



1994 



2010 



Liang and Qian, 2010; Cao et al. 2010 Zhou and Qian 2011). It has been constantly 



asked whether a general formulation exists. 

In the present work, we aim to provide a general and self-consistent description for 
the biological evolution, by considering it as a stochastic process and quantifying the 
evolutionary landscape as a potential function. We begin with one of the most classical 
models in population genetics, the Wright-Fisher model. Our construction comes prin- 



cipally from the novel decomposition of system's dynamical components (Kwon et al. 



2005), and coincides the Boltzmann-Gibbs distribution if the equilibrium distribution 



exists. It aims to capture system's middle- and long- term dynamical features. It treats 
all dynamical factors with unified viewpoints and approaches. The same ideology has 
already been used in studying phage-lambda genetic switch (Zhu et al. , 2004), Darwinian 



evolution (Ao, 2005) and cancer (Ao et al. , 2008). 



When studying the population dynamics with the constructed landscape, we en- 
counter a problem of estimating the waiting time for a stochastic escape from an infinite 
potential peak. For systems with multi-stable states, the issue of calculating the expected 
time/rate to transit between alternative peaks (or stable states) on a landscape has drawn 
lots of attentions in physics (escape rates), chemistry (reaction rates), and evolution (sub- 
stitution rates). A widely used approach to calculate the escape time is derived from the 
diffusion process, first introduced by Kramers (1940). It applies when the equilibrium 
distribution of probability density, near one or other peak, is approximately Gaussian. 
The escape time was also calculated using the backward diffusion equation in population 
genetics and other areas (Gardiner, 1985 Ewens, 2004 van Kampen, 2007) under similar 
approximations. The results estimate an exponential dependency of the escape time on 
the valley depth (barrier height in physics) on the landscape, called a "Arrehnius factor" . 
In the present model, however, this classical results cannot be applied as the potential 
landscape would generate infinite-high peaks under weak mutation. The Gaussian ap- 
proximation is invalid as the local maximum on the equilibrium distribution diverges to 
infinity. In this article, we develop a new way to analytically approximate the results 
around the divergent states, which provides a complete answer for the bi-stable dynamics 
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of the present model. 

The second approach for estimating the transition rate in population genetics is by 



calculating the probability of fixation of a single mutant in a finite population (Kimura 
1957, 1962). The rate of substitution (of a new mutation), then, is obtained as the 



multiplication of this probability and the average number of mutations that enter the 
population each generation (Gillespie 1998). This approach is valid under weak forward 
mutation, but is not generally applicable for more complex dynamics (e.g. when the 
backward mutation is not negligible). On the other hand, it provides a rather simple 
expression for the substitution rate, but the result is precise only under certain limits 
(e.g. when mutation is very weak in the neutral case). 

A third approach proposed by Barton and Rouhani (1987) is to calculate the leading 
(non-zero) eigenvalue of the diffusion equation, which gives the flux between different 
equilibria. Their method can be used when the forward mutation is weak, and allows the 
calculation under the joint effects of other factors (like backward mutation). However, 
their method fails when the selection is so weak that there only exists one peak on their 
"deterministic equilibrium" . The deterministic equilibrium fails to show the (long-term) 
bi-stability of the system in certain cases as the random genetic drift was not considered. 
Our U-shaped potential landscape, instead, is found better capturing system's dynamical 
essence and our results can be applied there. Also, their solution can be analytically 
approximated only under certain limiting conditions (very strong and very weak muta- 
tions); in the middle cases between the two limits, they need numerical computation 
for the integral term. Our results, instead, give analytical approximations for the whole 
parameter ranges that maintain the system's bi-stability. 

Other efforts have also been put in making analogies between statistical physics and 
population evolutions. Iwasa (1988) introduced a measure of statistical entropy and de- 
fined a free fitness function. This fitness function monotonically increases with time and 
takes its maximum at equilibrium. Barton and Coe ( 2009[ ) pushed this further and de- 
fined three analogous concepts from thermodynamics: a (negative) energy for selection, 
a potential function for mutation, and an entropy for random drift. Our present work, 
instead, integrates all above factors into an evolutionary potential landscape, and treat 
them with unified viewpoint and approaches. Other factors (e.g. frequency-dependent 
selection) can be readily added into our formulas. Barton and de Vladar (2009) used the 
analogy to thermodynamics to study the evolution of macroscopic quantitative traits, in 
a way that is independent of the microscopic details. However, there are cases (e.g. when 
mutation rate is low) that their maximum entropy approximation would not work. Also, 
their method requires a normalizable equilibrium distribution. They used the normaliza- 
tion constant as the generating function for macroscopic variables, which plays a "major 
role" in their calculations. We will show that our method is valid in the whole parameter 
space (even when weak mutation makes the potential function diverge), and does not 
require a normalizable equilibrium distribution. 

The present article is organized as follows. In Section 2, we analytically construct 
the potential landscape and specify it in the Wright-Fisher model; we analyze its basic 
structure and classify the mutation-drift models into several types; we discuss the uphill 
and downhill movements on the landscape, and the associated timescales. In Section 3, we 
calculate the first transition time between any two states and use it to study the average 
time to escape from an infinite potential peak, with the enlightenment of the landscape 
configuration. In Section 4, we further discuss the escape-time calculation and compare 
our results with previous efforts. We then show how further biological insights can be 
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gained from our framework by discussing three related issues: mathematical condition for 
complete fixation and absorbing boundary, mismatch of directed forces, and a comment 
on the misinterpretation of the stochastic tunneling effect. 



2 Model, Construction and Analysis 



2.1 Landscape construction 

The 1-D Wright-Fisher model considers the evolution of a diploid population at one locus. 
Assume that the generations are non-overlapping and the population size N is big enough 
for the continuous diffusion approximation. Denote the interested pair of alleles as Ai and 
A2 with respective frequencies x and 1 — x. Let p{x, t) be the probability distribution that 
Ai frequency is x at time t. The (1-D) diffusion equation for the continuous Wright-Fisher 



model is given by (Kimura, 1964 Ewens, 2004) 



dtp{x,t) = -dl V{x)p{x,t) 



- dx 



M{x)p{x,t) 



where M{x) is the average change of Ai frequency per generation, corresponding to cer- 
tain deterministic factors (e.g. mutation, migration and selection). V{x) is the variance, 
corresponding to random drift or other types of noise (e.g. fluctuated selection). 

In one of the present authors' work (Ao, 2005), a desired form of Eq.([T]) is expressed 
as follows: 

dtp{x,t) = dx eD{x)dx - f{x) p{x,t) . (2) 
The equivalence to Eq.([T]) is established by the following parameter transformations, 

fix) = M{x) - V'{x)/2 , (3) 
eD{x) = V{x)/2 . (4) 

where / is the directed force (we compare / and M in Section 4.4), D the diffusion term, 
and e the diffusion intensity. We take e = 1 in the present case. The motivation to take 
the form in Eq.(|2]) is simple: in cases where the detailed balance condition is satisfied, 
the potential landscape $ can be directly read from it (Ao, 2005). In 1-D case, it is 



^(x) 



fjy) 



dy 



(5) 



Here the sign is taken opposite to that in physics. 

On the other hand, by defining the probability current 

J(x) = f{x)p{x,t) - eD{x)dxp{x,t) , 

the equihbrium distribution for Eq.([2]) can be obtained under the natural boundary con- 
dition (Gardiner, 1985) 

0, (6) 



and the solution is 



p{x,t 



00 



J{x = 0,1; 
1 

— exp 



m 

eD{y) 



dy 



(7) 
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where Z is the normahzation constant 



p(x, t = oo)dx . 

Combining Eqs.(|5| [7]), there is the Boltzmann-Gibbs distribution (if Z < oo exists), 

p(x, t = oo) oc exp [$(x)/e] . (8) 

To specify Eq.([5]) in the Wright-Fisher model, we consider the joint effects of genetic 
drift, mutation, migration, and selection ( Felsenstein , 1976; Gillespie, 1998 Ewens, 2004). 



The average change of Ai frequency satisfies 

M{x) = —fix + z/(l — x) + m{Q 



x) + 



x{l — x) du 
2uJ dx 



(9) 



where p is the mutation rate from Ai alleles to A2 alleles per generation, and u is that 
from A2 to y4i; m is the fraction of immigrants from neighboring groups each generation 
with Ai frequency Q; cJ is the average fitness of the current population (a function of x). 
Under random genetic drift, 

V{x) = x{l- x)/2N . (10) 

2N is the number of alleles at the interested locus in the population gene pool. The 
potential landscape defined in Eq.(|5]) is then (consider the transformations in Eqs.([3| |4]) 

and set e = 1) 



$(x) 



lnx(l-x)+4A^ 



z/lnx+/iln(l- 



-x] 



+ 4:N 



niQ \nx+m{l—Q) ln(l— x) 



+2A^lncJ. 
(11) 

The four terms correspond respectively to the contributions of drift, mutation, migration. 



and selection. It is worth noting that $ is not proportional to 2N in Eq.(ll), so 1/2N 
may not be simply equivalent to the temperature in thermodynamics, which should be 



separated from the potential function in a Boltzmann-Gibbs distribution (see Barton and 



Coe (2009)'s discussion). 



We note that Eq.(^lip gives a general description for the system's dynamical elements. 



instead of defining different scalar measures for different biological factors as in |Barton and 
Coe (2009 )'s work. It provides unified descriptions and treatments for the system, and the 
associated approaches can be easily extended for other forms of biological factors. Second, 
Eq.([5]) is a dynamical construction, coming directly from the dynamical components of 
the system. Its validation does not depend on the normalizability of the equilibrium 
distribution (that Z is finite; see discussions of the pure drift case in Sections 2.2 and 4.2), 
which plays a "major role" in certain mathematical approaches (Barton and de Vladar 



2009). Its dynamical essence is on the ratio of potential values between different states. 
Nevertheless, when the equilibrium distribution is normalizable, it coincides with the 
definition of (negative) potential energy in a Boltzmann-Gibbs distribution in Eq.([8]). 
Third, our formulation is valid in the whole parameter space. We will show in Section 3 
that even divergent peaks on the equilibrium distribution can be handled in our framework 



(under weak mutations, compared to Barton and de Vladar (2009)'s work) 



2.2 General analysis 

We are interested in the evolutionary states which correspond to the peaks and valleys on 
the constructed landscape. The peak/valley here denotes the local maximum/minimum 
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on the landscape. It does not necessarily correspond to a most favorable state for selec- 
tion, but represents a (locally) most probable state as a result of the combined effects 
of all evolutionary factors. Likely, the uphill movement toward a local maximum on the 
potential landscape is not necessarily an "adaptive movement", and selection may be 
irrelevant. 

The first type of local maxima/minima is the system fixed points. If defined as x = a, 
it can be obtained by calculating 

= a) = f{a)/D{a) = . 

So x = a is the zero point of /. It is also a fixed point on the equilibrium distribution by 
Eq.([8]), if Z < cxD. Its stability (corresponds to a peak or valley state on the landscape) 
is determined by the sign of 

^"U = a) = - f(^)D'ix) 



The second type is the boundary states, namely $(0) and $(1), which may also be the 
local maxima/minima. An interesting thing in the Wright-Fisher model, is that genetic 
drift results in singularities at the boundaries. We show this in a simple system under 
mutation and drift. Specify Eq.([3]) with mutation and drift in Eqs.([9, 10), 

f{x) = -^x + z/(l - x) - (1 - 2x)/4:N , 
= F{x — a) , 

where we denote 

F = (1 - 2Nfi - 2iVz/) /2N , (12) 
a = {l-ANiy)/{2-ANn-ANu). (13) 

F measures the scale of the directed force /(x); a is the zero point of f{x) and the fixed 
point of $(x). The diffusion term D(x) is, by Eqs.(|4, 10) (setting e = 1) 

D{x) =x{l-x)/4:N . (14) 

The evolutionary potential landscape is then, by Eq.([5]), 

$(x) = (4A^z/ - 1) In X + {4:Nfi - 1) ln(l - x) . (15) 

We have $(0) = $(1) = oo under 4Ni>,4N^ < 1. The equilibrium distribution is (if 

Z <oo) 

p{x,t = oo) oc x^^''-^(l - x)^^^-i . 

Two divergent peaks emerge at x = 0, 1 for ANv^ANp < 1. Nevertheless, it is normaliz- 
able when u, fi > 0: 

= B{ANiy,ANij,) . 

The normalization constant has the form of the Beta function. When /i = or z/ = 0, it 
approaches infinity. = z/ = corresponds to the pure drift case, in which a population 
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will be fixed at either homozygotic states (x = 0, 1). The final distribution would depend 
on the initial probability distribution, and their is no equilibrium distribution (see also 



Section 4.3). Nevertheless, even under this setting, the landscape construction in Eq.(15) 
is still valid and changes continuously to be symmetrically U-shaped as /i, z/ — )■ 0. 

We may further define a partial normalization constant for each attractive basin, if 
< a < 1: 



a 



Zi= / x^^'^'^(l-x)^^'^-Mx , (16) 
Jo 

Z2= [\''"'-\l-xY^>'-'dx . (17) 



Each of them gives the probability ratio for observing a population in that attractive 
basin at equilibrium. Obviously Zi + Z2 = BlANujANfi); Z\ and Z2 are known as the 
incomplete Beta function. 



By calculating the fixed point and boundary values for Eq.(15), we classify and sum- 



marize the potential landscapes under mutation and drift in Table 1. 



Table 1: Classification of the mutation- drift model with changing parameters. "Mut- 
rate" means rate of mutation in certain directions (yU for Ai A2, for A2 — )■ Ai). Here 
we consider a fixed population size A^, so the system only changes with u, fi. x = a is 
the unique fixed point in the system (as / is linear) defined in Eq.(13). It's stability is 
determined by the sign of Eq.(12). The boundary values of $ is directly calculated from 
Eq.(15). The landscape configurations are determined by above factors. 



Mut-rate /i 


Mut-rate u 


Fixed point a 


Stability 


$(0),$(1) 


Landscape shape 


(0,l/4iV) 


(0,l/4iV) 


(0,1) 


Unstable 


+00, +00 


U-shaped 


(l/4iV,l) 


(l/4iV,l) 


(0,1) 


Stable 


—00,— 00 


Reverse U-shaped 


=1/4N 


=1/4N 


N/A 


N/A 


0,0 


Flat 


(0, l/4iV) 


(l/4iV, 1) 


Outside (0,1) 


N/A 


-l-oo,— 00 


Left-skewed 


(l/4iV,l) 


(0, l/4iV) 


Outside (0,1) 


N/A 


— oo,-|-oo 


Right- skewed 



2.3 Evolutionary timescales and bi-stable systems 

Evolutionary processes, visualized as movements on the potential landscape, can be clas- 



sified into two fundamentally different types: uphill and downhill processes (Zhou and 
Qian 2011). The uphill evolution is mainly driven by the directed forces. The downhill 
dynamics is often considered of stochastic essence (accumulations of rare events). This 
classification immediately suggests the separation of different evolutionary timescales. 
For example, in a bi-stable system, the first timescale Ti describes how long it would 
take to approach the local potential peak by the directed forces (may be interpreted as 
system's middle-term behaviors). It can be estimated as 

Ti~l/|/| , 

To specify Ti in the Wright-Fisher model, we choose a simple form of selection: Let s be 
the selective advantage of Ai over A2 (s <C 1), such that the average rate of change in x 



per generation (by selection alone) is (Kimura, 1964) 

M, = sx(l~x) . 
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The directed force then takes the form 



/(x) = sx{l -x) - fix + u{l - x) - {I- 2x)/AN . 



Note that in Eq.Q migration is mathematically equivalent to mutation by taking z/ = 
mQ,^ = m{l — Q), so we only discuss mutation here and for the rest of the work. The 
first timescale is (approximate sx{l — x) by sx/2 on average) 



2N 



\l-2N{n + u)±Ns 



-0{1) . 



(19) 



When 1/2N ^ s, (fi + u), the timescale is 2NO{l) (a well-known result for the pure drift 
case), and the evolution is dominated by drift; when (yU + z/) ^ l/2N,s, the timescale is 
(/i + z/)~^0(l), and the evolution is dominated by strong mutation (not likely to happen 
in natural systems); when s ^ l/2iV, [fi + z/), the timescale is s~^0{l) (a characteristic 
time for the selective sweep), and the evolution is dominated by selectional effects. These 
interpretations remind us of the characteristic times for certain evolutionary forces to 



1998). 



have significant effects on the genetic structure of the population (Gillespie, 

Note that Ti is salient when |/| > \D\, as the directed movement dominates the 
stochastic effect (neutral diffusion); when / ~ 0, the stochastic effect will dominate, and 
Ti becomes meaningless; when < |/| < \D\, the stochastic and directed factors are 
comparable. This corresponds to the change of skewness of the landscape configuration; 
note that f/D = $'(x) measures the slope of $. Also, Ti is a rough characteristic 
timescale for the uphill movement. The evolution rates should be different at different 
states. Near the fixed point, f{x) becomes comparable or even smaller than eD{x). The 
"last step" towards the fixed point may be manifested by the stochastic effects. 

The second timescale T2 is the time takes to reach the global potential peak (may be 
interpreted as system's long-term behaviors), and is dominated by the waiting time of 
escaping from a potential peak. The downhill movements in the peak-peak transition are 



often considered driven by stochastic factors (though see Waxman and Gavrilets (2005)'s 



work). Many efforts in literature tried to calculate the expected time/rate to escape from 
an attractive basin in specific contexts (Kramers 1940 Gardiner} 1985 [Lande 1985 



Barton and Rouhani| 1987[ Kriiger, 2010). A classical formula for the escape time [so 

T2 ~ Ti exp (A<l>) . (20) 



approximately T2) is 



Here A$ = ^{peak) — ^{valley) , is the valley depth on a landscape. The exponential term 
exp (A$) is often called the "Arrhenius factor" (Kramers, 1940); it shows the separation 



of the two timescales. However, in Eq.(15) and Table 1, there exist two infinite potential 
peaks at the boundaries (x = 0,1) under ANv^ANjji < 1 (which is also illustrated in 
Figure 1). As $(2; = a) is finite, we have 

A$ = $(0) - $(a) = 00 . 



Then by Eq.(20), the average time to escape from the potential peak at x = becomes 



To 



00 



The infinity may not be a good estimation for the escape time in all such cases with 
infinite potential peaks. Mathematically, under ANv = AN ^ = 1, the potential landscape 
becomes fiat (A$ = 0) and we obtain T2/T1 ~ 1; for any AN AN 11 < 1, however, we 
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obtain T2/T1 ~ 00. The escape time does not change continuously with ANu and 4N^ 
near 1, which is unexpected. Biologically, we expect that the evolutionary state would 
eventually reach an equilibrium distribution (though this process may be extremely slow), 
and for any finite-size system peak of the distribution would probably spread around. In 
the next section we show how we expand the classical formula and obtain analytical 
results for the escape time from infinite potential peak. 



3 Escape Time in Infinite Potential 

Without loss of generality, we only discuss the transitions in the direction — )■ 1 in this 
section. We start with the case of mutation and drift. First, we study the first transition 
time between any two states xq (near 0) and xi (near 1). We expand the integral terms in 
the classical formula near the divergent state x = and get an analytical approximation 
for the result. We show how the average time to escape from x = can be derived from 
this first transition time. We discuss the convergence of our expansion under various 
parameter settings. We then consider the effects of selection on the escape time by 
discussing two examples. 



3.1 First passage time 

We start with the mutation-drift case. The potential landscape is given by Eq.([l5|). To 
ensure a bi-stable system, we impose ANu, AN^ < 1; from Table 1, there is $(0) = $(1) = 



+00. The Arrhenius factor in Eq.(20) would give infinite time estimation for escaping 
from Xo = 0, and the solution would change very quickly with xq near 0. The question 
here is: Does the escape time really diverge with the infinite potential? Equivalently, 
does the infinity of equilibrium distribution at certain states imply the impossibility of 
peak-shifting events? If not, how does it depend on the parameters of the system? 

As the Gaussian approximation fails, we come back to the standard formula for the 



first passage time, derived from the backward diffusion equation (Gardiner, 1985). Our 
first consideration of the escape time is enlightened by the landscape configuration in 
Figure lb, where a potential valley x = a separates two attractive basins. It implies 
different system behaviors in different frequency intervals (0, a) and (a, 1). A large portion 
of the escape time is often shown to be spent in the small region near the valley (or barrier) 



state (Kramers, 1940). 

We study the first passage event through the valley point a; = a to some state xi > a, 
starting from some state xq ~ in the attractive basin (0,a). The average first passage 
time T{xo) satisfies 

[fix) + eD'{x)]d,T{xo) + eD{x)dlT{xo) = -1 . 

To study the time to exit (0,xi) from xi, we set the boundary conditions as (x = 
reflecting and x = Xi absorbing) 

d,T{x = 0) = , 
T(x = Xl) = . 



The solution is 



T(xo^Xi) = y ■^^^exp[-<l>{y)]dy J exp [(l>{z)]dz . (21) 
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Here $ is just our evolutionary potential landscape in Eq.(15). Previous approximation 



methods that generate the Arrehnius factor is mainly established on the following two 
assumptions (Kramers, 1940 Gardiner, 1985): (1) The main contribution to the integral 
in Eq.([2l|) is due to a small region near x = a (so the main time of escape is spent near the 
region). (2) Assume Gaussian-like equilibrium distribution around x = (so $'(x = 0) 
equals zero) . In the present illustrated by Figure la, the equilibrium distribution 

is not typically sharply peaked at x = a; there may be a "fat" potential valley rather 
than a "thin" one. The bulk of the integral may span over a big region around x = a, 
and thus assumption (1) may not be true. Biologically, a population may wonder around 
the valley state for a considerable amount of time instead of going straightforwardly to 
the other potential peak. More importantly, the assumption of Gaussian-like equilibrium 
distribution fails near x = 0: Not only is <l>'(x = 0) non-zero, it actually approaches 



tuting Eqs.(14, 15), and there is 



infinity. The classical result Eq.(20) cannot be applied for this situation. 



Here we specify Eq.(21) under weak mutation {ANiy,ANfi < 1) and drift by substi- 



T{xo ^ xi) = AN r y-^^'^il - yY^^'^dy F z^''^-\l - z^^^-'dz . 

JxQ Jo 



(22) 



The second integral is just the integral of p{x,t = oo) in [0,y], having the form of 
incomplete Beta function. No previous expansion was taken on this divergent equilibrium 



distribution, though Barton and Rouhani (1987) claimed that their eigenvalue method 



can be applied in the Gamma-like divergent equilibrium (see Section 4.2). We describe 



our expansion of Eq.(22) as follows. 



Under < z < y < Xi, there isO<l — a;i<l— 1/<1 — 2;<1; the expansions on the 
exponential terms of 1 — y and 1 — z would converge near 0. First we expand (1 — 2;)^^^"^ 



with respect to z in the second integral of of Eq.(22) (denote it as I{y)) 



y 



y 



z 



4Nu 



-Eri(^)- 

n=l k=l ^ ^ 

4iVz/^J-J-l k Jn + ANu 

n=l k=l ^ ^ 



dz 



(23) 



The convergence is obvious given < y < Xi < 1. Substitute and expend (1 — y) ^^'^ 
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in the first integral of T(xo — )■ Xi), we obtain 



T(xo — > Xi] 







A AT 






xo 












xo 


Xi - 


Xo 


V 




X\ — 


Xo 



oo n 

n=l k=l 

_ oo n 



A; 



1 ^AA--4iVfi\ 



1/ 



n=l fc=l 



A;-l + 4A^/i\x^ 



n+l 



A; 



4iVz/ ■ ^ J-J- V A; 
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(24) 



With Xl, Xo < 1, the convergence of this expansion is obvious. We discuss its convergence 
under two hmiting conditions: 

(1) Under x — )• and — )■ 0. For any non-zero z/, our result for escape time in Eq.(28 ) 
is not sensitive to xq near 0, and we can take xq — )• 0. It is of the essence that the 
incomplete Beta function in Eq.(23) is of near-linear dependence on the integral interval 
near the boundary x = 0, given u > 0. For z/ = 0, however, the expansion in Eq.(23) will 
not be valid. We have instead 



liy) 



n=l k=l ^ ^ 



The leading term changes from a polynomial of order ~ y^^ to the logarithmic scale 
(~ Iny). The result becomes sensitive to the value of xq near and approaches infinity 
at Xq = 0. To ensure the convergence of the escape time as Xq — )■ under weak mutation, 
we should first take Xq — )■ 0, and then we can set z/ — 0. Nevertheless, as discussed above, 
the result T(xo — ?■ a) — )■ oo is consistent for any way of taking z/ — )■ 0. 

(2) Under x — )■ 1 and /i — )■ 1/4A^. When 4A^/i = 1, the expansion of (1 — y)~^^^ in 
Eq.(24) would not converge for x\ = 1, as the series 

n=2 k=2 ^ 



k 



n 



(25) 



would then become a harmonic series (with logarithmic divergence rate) in maths. We 

yy^dy 



return to Eq.([22j) for iNfi = 1: 
T(xo — )■ Xl) = 



4N 



Xl 



dz 



V 



-Mn 



Xo 



Xl 



As Xl — )■ 1, this result approaches infinity. This is also illustrated by the vanishing bi- 
stability of the system (see Figure 3 that x = 1 becomes then a potential valley; also see 
Section 3.2 for more details). To ensure the validation of the expansion, we require at 
least that Xi < 1 or 4A^/i < 1. 
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3.2 Escape time 



Biologically, we want to know the average time to escape from x = 0, after which the 
population is not expected to return immediately. The escape time can be derived from 
above first passage time by setting xq ~ 0,a;i ~ 1 (see also Section 4.1). This statement 
comes from our following observations on the landscape configuration shown in Figure lb 
{4:Nh',4:Nfi < 1): Near x = a, the landscape valley is flatly distributed and a population 
may wander there for some time before approaching either boundary state; and for any 
a < Xi < 1 (xi not near 1), it may still come back to x = with considerable possibilities 
after reaching xi; near x = 1, the landscape peak grows sharply, so that a population 
would be caught stable there (not easily returns to x = 0). Thus the escape time can 
be approximately considered as the average time taken to reach another stable state 



(substitution time in population genetics). We take xq — )■ 0, xi — )■ 1 in Eq.(24): 



r(0 ^ 1) = 4iV y-^'^'^il - vY'^^dy f 
Jo Jo 



(26) 



1 1 

- + - 



4iV 



... A A - l + 4Ar/i\ 1 

2iv/i+4iv/.^n — I — -) 

n=2 k=2 ^ ^ 



n + 1 



2(1 + 4A^z/) 



f:n(^) 

n=2 k=2 ^ ^ 



[n + l)(n + ANu) 



(27) 



In Eq.(24), as Xi — Xq increases, the convergence of the geometric series becomes slower, 
and the leading term approximation becomes less accurate. If we set xq = 0,xi = 1, 



the sum of the geometric series Yl^i diverges. The expansion in Eq.(27) does not 



necessarily diverge, though, because the geometric terms are multiplied by other series 



factors. For example, the second integral of Eq.(27) converges with rate 1/n^. The 
first integral needs more detailed discussions on the value of AN ft. It has two effects 



on the expansion: as a series factor it determines whether Eq.(25) converges and as a 
pre-coefficient it gives the weight of this series (if it converges). For the first effect: When 
ANfi > 1, Eq.(25) obviously diverges with the rate of at least l/n; when 4Nfi < 1, we 



prove in the Appendix B that it converges. For the second effect of 4A^/i as a pre-factor of 
the first series (under ANfi < 1), it gives the weight of the remainder terms. If ANfi <^ 1, 
the remainder terms are negligible and the leading term approximation is accurate; if 
4N^ ^ 1, the remainder terms become important and the first several terms in Eq.(27) 
cannot approximate the escape time very well. 



We make a further comment on Eq.(27). Under ANu <^ 1, the scale oi l/u would be 
much bigger than the remaining terms (of order 2N). The interesting thing is that even 
though the expansion is taken on a, it naturally separates the two distinct timescales. 



and provides a good analytical approximation for the escape time in Eq.(22). If instead 



4A'" and ly ^ are comparable (under the constraint ANu < 1), the expansion in Eq.(27) 



is still valid, but there would be no clear separation between the leading term and the 
rest. It implies a mixed timescale of Ti and T2, illustrated by the near-fiat slopes on the 
constructed potential landscape (Figure 3a is an example for comparable fi~^ and 4A^). 



Assume ANv^AN^ <^ 1 and take the leading terms in Eq.(27), we get 



T(0 ^ 1) 



1 + 2A^/i 



+ 2A^C(1) . 



(28) 
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Comparisons between the numerical solution of the classical formula Eq.(22), our ana- 



lytical approximation in Eq.(28), and the results in the discrete model (see Appendix A) 



Eq.(Eq:Al-l) is given in Figure 2. 



From Eq.(28), the escape time increases with and /x; this is consistent with the 
biological interpretations. For increasing fi, the stronger backward mutation makes the 
forward transition more difficult; for 4Nfi <^ 1, the effect is near-linear; for ANfi > 1, it 
makes the transition to x = 1 impossible, as discussed in Section 3.1. This is directly 
illustrated by the change of landscape configuration in Figure 3. The fixation of Ai 
becomes phenomenally impossible because the attractive basin (a, 1) vanishes and x = 1 



becomes an unstable state (potential valley) then. By Eq.(22) we get 

,1 



4N 



XI 



y 



dy 



xo 



1 



-dz 



y 1 

This form suggests an exponential dependency of the transition rate in the direction 
Xo — )■ Xi on AN IX. For increasing A^, it has effects from two sides: increase the number 



of new mutants each generation and decrease the intensity of random drift. By Eq.(28), 
the net effect is to increase the escape time, though in a minor order if ANv ^ 1. 
If further taking /i 



0, Eq.(28) becomes 
T(0 ^ 1) 



V 



(29) 



The result is approximately independent of the population size A^. It reminds us of the 
well-known formula for the rate of substitution for the neutral mutants: 2Nv x 1/2N = v 



Gillespie , 


1998 


); Eq.( 


29 



characteristic timescale for mutation, but coincides with it under the limit 4A^i/ ^ 1. 
The driving force for the process is mainly random drift (1/2A^ 3> v\ without which the 
stochastic escape cannot happen, even though the effect of A^ cancels out under the limit. 



Like in Eq.(24), the escape time in the direction 1 — )■ can be estimated as 

vY^'^^dy 



T(l -> 0) = 4A^ 

= 4A^ 
1 




1 

- + - 
II /i 



yY^'^'dy 

oo n 



z 



\4:Nu-l 



n=2 k=2 



k - l + ANiy\ 
k ) 



n 



dz 



+ 



1 - ANv 
2{l + AN^i) 



n=2 k=2 



ANp\ 1 

) {n + l){n + AN^) 



k 



Transitions in the two directions will eventually balance each other, and the global equi- 
librium distribution will be established. 



3.3 Effect of weak selection 

The general formula for the first transition time when there is mutation, drift and selection 
can be obtained by substituting Eq.(ll) into Eq.([2l|) (we neglect migration because it is 
mathematically equivalent to mutation) 

pxi py 
T(xo ^X,)=m (1 - y)-^^>^y-^^-^ Piy)] dy / (1 - z^^^^-^ z^^-^'^ [cJ(^)] dz . 
Jxo Jo 
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If we can expand [aJ(?/)]~^^ and [cJ(2;)]^^ near 0, an analytical approximation can be 
obtained by combining it with the previous expansion results described in Section 3.1. 
For example, suppose the simple form of selection defined in Eq.([l8)) (s is the selective 
advantage of Ai over s -C 1); the evolutionary potential landscape becomes 



= {ANu - 1) Inx + {4Nfi - 1) ln(l - x) + ANsx 



(30) 



The valley state a is defined by taking = a) = 0. For maintaining a bi-stable system, 
we set 1/4N > /x, (see Figure 4a). The first transition time is obtained as 



T(xo xi] 



/■Xl 

AN / e" 

J XO 



4Nsz 



(1 _ ^)4^^-1^4iV-lt/^ . 



To expand e~^^*^ and e^^*^ near 0, we assume ANs < 1. Combine the expanded polyno- 
mial series of the exponential terms with the previous results, we obtain the escape time 
(set = 0, xi = 1 and assume ANu ^ 1) 



T(0^1)^4A^y" 1+ (^4:Nfi-4:Ns^y 
1 + 2iV(^ - s) 



1 l-4Nu + 4:Ns 

+ — ^^TTT y 



ANv 



dy 



(31) 



The convergence is obvious given 1/4A^ > /x, z/, s. From this expression, the selective 
advantage decreases the escape time on a near-linear scale. The result is consistent to 



the rate of substitution k (Gillespie, 1998) under the same parameter settings (/x = 0, s <^ 
l,ANs < 1): 



k 



1 - 

71^ X 2A^z/ , 



2s 



4A^s - 8{Ns)^ 

V 



X 2Nv 



1 - 2Ns 



just the reverse of Eq.(31) if we take yU 
under 4A^s < 1. For 2A^s > 1 



0. The approximation in Eq.(31) is accurate 
it obviously fails as Eq.(31) would give T(0 — )■ 1) 



< 0. 

Another comment for Eq.(31) is that the escape time is more sensitive to the mutation 
rates than the selection intensity if both effects are weak [ANv^ANs < 1). When u 
changes from 1/4A^ to 0, T2 increases from 1 to 00; when s changes from 1/4A^ to 0, T2 
is just changed by a factor of 1/2. It is consistent with the observation that "the low 



rate of mutation limits the ability of a species to reach the most fit sequence" ( Gillespie 
1984) for this setting. 

Another example of selection is to take Ms{x) = — sx(l — a;)(l — 2x) and /i = i^, where 



s is the fitness deficit of the heterozygote relative to the homozygotes (Lande, 1979). We 
have 



$(x) = {ANfi - l)lnx(l 



x] 



ANsx + ANsx^ 



(32) 



which maintains two "deterministic equilibria" for s > 4/x (Barton and Rouhani, 1987). 



Barton and Rouhani (1987)'s method to calculate the peak-shifting rates failed for this 



model under s < 4/x, as their "deterministic equilibrium" becomes uni-peaked then. We 
show in Figure 4b that the system is still bi-stable under this parameter setting, and our 
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expansion can treat this case. By taking Eq.(32) into Eq.(21), the transition time from 
to 1 is obtained as (assume 4Nfi <^ 1, < 4A^s < 1 and expand the exponential terms 
hke in the previous example) 

T(0 -> 1) 

= AN f e^'^'y-^^'^y^il - y)-^''^y-^''^'dy T e-4A^-+4Af-^(i _ ^)4^a^-i^47v^-Mz , (33) 
Jo Jo 

^AN f (1 + ANsy - ANsy'^){l + ANfiy)y-^^^ f\l - ANsz + AN sz'^)z^^^-^dz , 
Jo Jo 

^ -(l + 2Nfi+-Ns) . (34) 
/i V 3 / 

The transition time shows linear dependency on the selection coefficient s, which is verified 
by the numerical solution in Figure 5. This is also consistent with the Monte Carlo 



experiment of Barton and Rouhani (1987) for ANs < 1. The expansion is exact for 



ANs <C 1, so is valid for s < fi. The error increases as 4A^s — )■ 1 because the high order 
terms of ANs become more important then. 

4 Discussion 

4.1 Escape time and first passage time 

Qualitative analysis of Figure lb in Section 3.1 (potential landscape with "fat" valley and 
sharp peaks) suggests a very different transition process from the classical escape problem. 



This is quantitatively verified by comparing the result in Eq.(24) to the classical form 



in Eq.(20), mainly the absence of the Arrehnius factor. The main contribution to the 
waiting time is not crossing the potential valley at x = a; instead, it shows a near-linear 
dependence on xi — Xq (under ANujANfj, <^ 1). The wanted escape time cannot be 
approximated by the first passage time through the valley point x = a to a nearby state. 

Let's do a thought experiment here. After passing through the valley point and 
appearing in (a, 1), a population may still wander around the valley state for some time, 
and may (even with considerable probability) come back into (0, a) and reach x = in 
Ti. Actually, if T2 ^ Ti, the escape process from the fixation of Ai can be equivalently 



considered as a series of Markovian passage events. In Eq.(15), under ANu, AN jj, <C 1, the 



landscape configuration is almost symmetric to a = 1/2. The first passage time through 



the valley point, by setting xq = 0,xi = a and taking AN 11, AN v — )• in Eq.(28), becomes 

T(0 ^ a) ^ {2u)-^ . 

The rates of uphill (~ 1/2N) and downhill (~ u) movements are well separated. Once a 
population reaches x = a in l/2i/, it has (on average) probability 1/2 to fall into either 
attractive basin and immediately (~ 2A^) reaches the corresponding potential peak; once 
falling back to x = 0, it will wait another l/2z/ to reach x = a again; it then again has 1/2 
chance to reach x = 1 or return to x = immediately. Assume this process continues. 
The expected time (also approximately the second timescale T2) for leaving x = under 
ANu, AN 11 <^ 1 can then be obtained as 

^112 1 n 1 

T ~ — x-H X — + ...H X h... , 

2u 2 2v 22 2v 2" 
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the same result as Eq.(29). 



4.2 Comparison with previous efforts 



Kramers ( |1940 ) derived the rate of escaping over a potential barrier from the diffusion 
equation, by calculating the probability current at the saddle point (our valley state). 
He assumed that the main contribution to the escape time is due to a small region 
near the saddle. His analytical approximation requires a finite barrier height, and that 
the equilibrium distribution near the potential well be approximately Gaussian. These 
assumptions generate an "Arrehnius exponential factor" as in Eq.(20) in the calculations. 
Other methods developed from the backward diffusion equation obtained similar results 
(Gardiner, 1985; Lande, 1985). In population genetics, however, the form of random 
genetic drift in Eq.(14) generally leads to non-Gaussian (even divergent) equilibrium 
distributions. The potential landscape is typically flat at the middle states and divergent 
at the boundaries. This constrains the application of the classical escape formulae in 
certain cases: The peak-shifting rates are often numerically approximated or can be 



analytically discussed only under very special settings in population genetics (Kimura 



and Ohta, 1969). Our expansion Eq.(24) taken on the escape formula Eq.(22) near the 



divergent state makes it possible to analytically approximate the results in such cases. 
It is proved to be valid for all the conditions where a bi-stable potential landscape is 
maintained. Our results show that the finiteness of the escape time does not require the 
finiteness of the potential peak (or barrier height), but corresponds to the normalizability 
of the probability distribution in that attractive basin (see Section 4.3). 

As mentioned previously, our result in Eq.(29) can be compared with that derived 
as the rate of substitution in population genetics (Kimura, 1957, 1962). First, the es- 



timation of neutral substitution rate v (Gillespie, 1998) only gives a characteristic rate 
for the fixation of Ax to take place. We cannot have a detailed knowledge about how 
the transition time depends on other factors of the system. An example is the effect of 
population size A^, though it may be of minor impacts under v <C 1/4A^. If instead v and 
1/4A^ become comparable, the simple estimation would not work any more. We've shown 
that our analytical approximation is valid for the whole range < 4A^z/ < 1; when 4A^z/ 
is not very small, the timescales of the uphill and downhill movements become mixed 
and the substitution time should be expressed by a more complex form. Second, the evo- 
lutionary transition time calculated in the present work can be more generally applied 
in population genetics than the substitution time, because all evolutionary factors are 
here treated in a unified and consistent manner. For example, it allows the existence of 
backward mutation (yU > 0), which may make the fixation probability of Ax incalculable. 
We've shown in Sections 3.1 and 3.2 that it has a big effect on the average time to transit 
between states. 



Barton and Rouhani (1987) proposed to calculate the leading (non-zero) eigenvalue of 
Eq. (l2j), which gives the flux between different equilibria. The authors claimed that their 
general expression of transition rate agrees with the two approaches (backward diffusion 
equation and rate of substitution) in the limits of very high and very low mutation, 
respectively. Our first comment is that, in this sense, our expansion also bridges the 
two approaches in the two extreme conditions — we make the first method applicable in 
the second situation. Second, their method failed under very weak selection s < 4fi. It 
requires the existence of two peaks on the "deterministic equilibrium", which does not 
necessarily capture the correct dynamical features of the system. For example, even in 
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the neutral case (s = 0, 4A^yU, 4A^z/ < 1), the equihbrium distribution (and our potential 
landscape) has two local maxima, implying the long-term bi-stability. We've showed in 
Section 3.3 that our result can be applied for the whole range of < 4Ns < 1. Third, 
part of their solution is numerical and can only be analytically expressed under limiting 
conditions [ANfi ^ 1 or ANfi ^ 1). Our method, in contrast, can be used for analytical 
results not only under extreme conditions, but in the cases with intermediate parameter 
values (when u and 1/4N are comparable). It naturally separates two distinct timescales 



from the complex form of Eq.(22), and gives a good analytical approximation for the 



escape time. Our result provides a complete answer for the present bi-stable models. 

4.3 Complete fixation, absorbing boundary and normalization 
constant 



By taking z/ = in Eq.(28), we see that the expected time to escape from x = becomes 

T(0 -> 1) = cx) . 

It implies that no escape can happen once trapped into the potential peak. If further 
assuming /i 7^ 0, we have T(l — )■ 0) < 00; there would be always probability flux in the 
direction 1 — )■ 0, unless the probability densities in (0, 1] all vanish. We can expect, after 
a (finite) long time (~ fi~^0{l)), all probability densities will be fixed at a; = 0. This is 
consistent with the equilibrium solution of Eq.([2]) (McKane and Waxman, 2007). Under 



/i > = 0, the peak at x = reduces to a Dirac delta function S{x)^ and the probabilities 
at other states all vanish; for z/ > /i = 0, this happens atx = l. li u = fi = (pure 
drift case), there will be two delta functions, whose relative weights depend on the initial 
distribution; the final distribution would be p{x, t = 00) = (1 — xo)S{x — 0) + xoS{x — 1), 
where Xq =< xp{x, t = 0) > is the initial average state, and there would be no unique 
equilibrium distribution. We call it "complete fixation" that a population starting from 
any initially state will finally be fixed at a specific homozygotic state (e.g. a; = 0). Such 
condition for the complete fixation can thus be derived from our escape time results. 



From the convergence analysis of Eq.(27), infinite escape time can also be achieved by 
setting /i > 1/4N] however, under this condition, the potential landscape becomes uni- 
peaked (Table 1 and Figure 3), and the escape issue becomes irrelevant. Similar is for 
setting N ^ 00, that all populations will deterministically move to a heterozygotic state 



a = u / {jjL + u) (Eq.(13)), instead of fixing at any homozygotic state. 

From above analysis, we find that x = becomes an absorbing boundary under z/ = 0. 
This observation is interesting because we initially set natural boundaries (Eq.do])) for the 



model, but later derive an absorbing boundary from the dynamical equation. There is 



no real conflict between the two settings. McKane and Waxman (2007) gave a general 



formula for the boundary condition in the diffusion equation (their Eq.(lO)), which can be 
readily obtained from our results for the escape time. We make still another connection 
to their work: The continuities of the diffusion equation (as a theoretical description for 
the system's dynamics) from (0,1) to [0,1] (including boundaries) and from z/, /i > 
to z/, yU > are verified by our demonstration of the convergence of the escape time. 
Our framework established on the evolutionary potential landscape, as discussed above, 
contains essentially the complete information of the frequency space and the parameter 
space in a consistent manner. 



In Eq.(22), the second integral I{y) is just the partial normalization constant in a 



region (0,?/). If we can approximately set y = a for J, the infinity of escape time would 
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directly come from the unnormalizability of p{x,t) in [0, a] {Zi = oo, defined in Eq.(16)). 
The condition for the escape time to diverge as x — )■ 0, is not the divergence of the equi- 
hbrium distribution (or the potential landscape) at that potential peak, but the infinity 
of the partial normalization constant in its neighboring region (attractive basin). The 
divergence of the probability density function is a necessary but not sufficient condition. 
We summarize above observations in Table 2. 



Table 2: Summarization of the relations among the (partial) normalization constants, 
the escape time, the fixation type, and the boundary condition. Zi and Z2 are the 
partial normalization terms defined in Eqs.(18, 19). "Fix on prob." means the fixation 
may happen at x = or x = 1 on probability (depending on the initial distribution). 
T(0 — 7- 1) and T(l — )■ 0) are the corresponding escape times. "Absorb-bound." means 
where the absorbing-boundary condition emerges. 





Z2 


T(0 ^ 1) 


T(l 


^0) 


Fixation type 


Absorb-bound. 


< 00 


< 00 


< 00 


< 


00 


N/A 


Neither 


= 00 


< 00 


= 00 


< 


00 


Complete (x = 0) 


X = 


< 00 


= 00 


< 00 




00 


Complete (x = 1) 


X = 1 


= 00 


= 00 


= 00 




00 


Fix on prob. 


X = 0,1 



4.4 Comparison between the directed forces: / and M 



One of the present authors (Ao et al. , 2007) proposed a new approach to treat the 



stochastic differential equation (SDE): 



dx/dt = /(x) + ^/2eD{x)^{t) 



(35) 



where ^(t) is the Gaussian white noise. In the one-dimensional case, this treatment (we 
may name it A- type according to Shi et al. ( 2011[ )) can be compared to those of Ito and 
Stratonovich by the following manner. If by Ito convention the diffusion term D{x) takes 
the value of x just before the jump, and by Stratonovich convention it takes half the sum 



of the values before and after the jump (van Kampen 1981 ), then by A-type convention, 
D{x) takes the value of x just after the jump. Eq.([3]) in fact reflects the relation between 
the directed forces of Ito's (M) and A-type (/) under the 1-D case. 

From Eq.^, M can be considered as the zero-noise limit of /; mismatch of the two 
implies the shift of system flxed point after introducing certain multiplicative noise (that 
D{x) is not constant). If D{x) is constant for x G [0,1], there will not be such extra 
term. The fact that multiplicative noise contributes to the directed force has already 



been known in the study of dynamical system (Ao, 2005). By taking the A-type integral 



the SDE Eq.(35) would generate the Fokker-Planck equation Eq.(|2|, and there will be no 
such mismatch issue. The mismatch of / and M vanishes as e — >■ 0; in the Wright-Fisher 
model, this is manifested by taking N ^ 00. 

By Eq.(|5]), / (instead of M) is directly related to the gradient of the potential and 
captures the property of the Boltzmann-Gibbs distribution (if the equilibrium distribution 
is normalizable). In the Wright-Fisher model, M(x) integrates all deterministic factors, 
such as mutation, migration and selection. It is deflned as the average increase rate of x 
(it can be obtained from Eq.([2]) that d < x > / dt =< M{x) >, see Appendix C). Under 
mutation and drift, M can be expressed as M(x) = — /xx + v{l — x) = Fm{x — aM), 
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where Fm = — (/i + i^), = z//(yU + z/) measure the scale and zero point of M. They can 
be compared to the similar form of f{x) in Eqs.(12, 13). By setting fi = u < 1/4A^, we 
obtain 



a = aM = 1/2 , 
F>0> Fm . 

Then a; = 1/2 becomes a stable zero point for M but an unstable point for /. In Table 
1, the landscape is U-shaped (so is the equilibrium distribution of p), and x = 1/2 is 
the global minima. It is consistent with the effect of / (at least qualitatively). If further 
taking fi = u = (pure drift), we obtain M = 0, F = 1/2N, a = 1/2. The first timescale 
is of order 1/F = 2N by Eq.(19), just the well-known result under pure drift (Kimura 



1964); the quantitative side is verified. 



4.5 A comment on the "stochastic tunnehng" 



Iwasa et al. (2004) studied a finite population of replicating cells of three types, with 
unidirectional mutations in two transition processes (type — )■ type 1 — )■ type 2). They 
termed a stochastic "tunneling" phenomenon that allows transition from the fixation of 
type cells directly to the fixation of type 2, without passing through the fixation of 
type 1. We make two comments here. The first is that their use of the term implies the 
existence of a potential barrier (or landscape valley). Conceptually, tunneling means a 
barrier-crossing process in a potential field. Second, their description of the phenomenon 
is misleading. "Tunneling" refers to the quantum mechanism which is classically impos- 
sible, but here the transitions are driven by the stochastic factors. The phenomenon 
appears because they tried to model the non-Markovian process (the fixation of type 
cells — )■ the fixation of type 1 — )■ the fixation of type 2) as a Markovian process. They 
failed when the timescales of the two transitions are comparable or the first transition is 
slower. It should be a process of climbing over a saddle point on the landscape, driven by 
the accumulative perturbations. The proper term here should be "thermal" activation. 
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Appendices 

A. Escape time in the discrete model 

To compare with the diffusion model and verify our Eq.([28|, we calculate the escape 
time in the discrete Wright-Fisher model. We consider a Master equation for the Moran 
model, which gives rise to the same mesoscopic equation Eq.Q as the Wright-Fisher 
model. The only difference between the two models is that the timescales differ (Blythe 



and McKane 2007). Here in our formulation, a generation in the Wright-Fisher model is 



2N times as long as that in the Moran model. 



Gardiner ( 1985 ) gives the solution of the mean passage time from x = to a; = 2A^ (x 



in this section denotes the number of Ai alleles in the population, instead of its frequency) 



2N-1 



{A.l] 



y=0 



=0 



X = 2N is set as the exit state. is the probability of transition x — )■ x + 1 per unit 
time; t~ is that for transition x — )■ x — 1. 0(x) is defined as 

When specified in the Moran model under mutation and drift, the parameterizations are 



(Blythe and McKane, 2007): 



r(n) 



z/ 1- 



n 
N 



N 



1 - 



N 
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Note that r(0 — )■ 2N) in Eq.(A.l) should be divided by 2N when making comparison 
with the Wright-Fisher modeL Simple analysis shows that u = makes 0(0) = 0, and 
then r(0 — ?■ 2N) = oo; this is consistent with Eq.(28). 



B. Convergence of Eq.(25) 



To study the convergence of the partial sum 

'k-l + 4:Nfi\ 1 



n=2 k=2 



k 



n + 1 



We use Raabe's test for series convergence from standard textbooks of real analysis. For 
< ANn < 1, we denote 



n 

k=2 



k-l+4:Nfx\ 1 



k 



n + 1 



Obviously c„ is positive for all n > 0. First, we have 

lim ^"''^^ = 1 . 



n-s>oo c„ 



{B.l] 



Then we calculate the Raabe terms 



J-, I Cn+i 
Rry = n\ 1 
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n+1 
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>- k=2 

n + 4A^/i 
n + 2 



k-l + AN II 1 
k n + 2 
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k=2 



k-l + AN fx 1 
k n + 1 



{ANfi - 2) 
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n + 2 



Here AN^ — 2 is a constant less than —1. By taking the limit n — )• oo, the Raabe term 
becomes 

lim i?„ = 4A^/i-2< -1 (B.2) 



under < ANfx < 1 (Thomson et al. 2008). 



The two conclusions in Eqs.(B.l, B.2) verify the convergence of the partial sum of S„ 



C. Movement in the first moment 

We study the transient movement of a population on the landscape. Assume a population 
has an initial state Xq. Then at time t = we have p(xo,to) = ^{^ ~ ^o)- In the first 
moment, the average speed of its movement is calculated as 



d < X > 
dt 



d 
dt 



xp{x, t)dx 



'x'^dx 
dt 
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The average < ... > is taken over the population's possible states; at the first moment it 
is the Dirac delta function centered at x = xq. By substituting Eq.Q, we get 



d < X > 
Jt 



dx 
' d 

X 



eD{x)^^^^~f{x)p{x,t) 







dx 



eDix 



dx 
dp{x,t) 



dx 



dx — 



dx 



dx 

f{x)p{x,t) 



dx 



.dp{x,t) ^ 



eDix 



dp{x,t) 
dx 



dx — f{x)p{x,t)x\l + / f{x)p{x,t)dx 



Note that the first and third terms in sum gives —J{x,t)\l, which is 0; the fourth term 
is actually < /(x) >; the second term is 



eD(x 



dp{x, t) 
dx 



dx = —eD{x)p{x,t)\l + / eD'{x)p{x,t)dx 

Jo 



Again the first term is as = 0, 1) = and the second is < eD{x) >. So in all 

d < X > 



dt 



--< fix) > + < eD'ix) > 
--< Mix) > 



This is just the definition of average transition rate in the diffusion equation. By imposing 
p(xo,to) = ~ 2^o)) we obtain the expected transient movement rate M(xo). 
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Figure 1: Equilibrium distribution and the corresponding adaptive landscape for a typical 
bi-stable system under mutation and drift. The parameter is set as 1/4N > u > fi (drift 
dominates mutation); the peak at x = 1 is higher than that at a; = on each figure, 
(a) As z/, /i — )■ 0, the probabilities would accumulate more sharply at the boundaries; 
when z/ = /i = 0, it develops two Dirac delta functions, (b) The U-shaped landscape 
shows a "fat" valley at the middle states and two sharp peaks at the boundaries. A 
population may wander some time at the middle (heterozygotic) states, but would be 
absorbed quickly into the boundary (homozygotic) state once appearing near one. As 
z/, — )■ 0, the landscape would continuously converge to a symmetric configuration; there 
would be no sudden change in the limit u = fj, = 0. 
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Figure 2: Comparisons among the numerical solution of Eq.(26) (solid line), analytical ap- 
proximation in Eq.(28) (dashed line), and results in the discrete model (crosses, Eq.( A.l[ )) 
in the range < ANu < 1. Other parameter settings are: /i = 0.00005, = 100. 
The leading terms in the expansion Eq.(28) approximates the exact solutions well for 
ANiy ^ 1; when ANu ^ 1, the two timescales are comparable and the remaining terms 
becomes important, so the error becomes major. The numerical solution of the dis- 
crete formula Eq.(A.l) matches well with the continuous results for the whole range 
< ANu < 1, though is found generally smaller than the latter. This may be explained 
by that a diffusion equation is the second order approximation of a master equation. 
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Allele frequency of 

Figure 3: Adaptive landscapes for strong backward mutations [ANu > 1). The attractive 
basin (a, 1) vanishes and the homozygotic state x = 1 becomes unstable, (a) ANfi = 1, 
the critical case; bi-stable system changes to be uni-peaked, and the valley point a locates 
at X = 1. The transition 0-^-1 becomes pure downhill; even if a population reaches x = 1 
with non-zero possibility, it would return quickly to x = 0. Escaping from x = becomes 
phenomenally impossible, (b) 4A^/i > 1, the valley point a locates outside (0,1), and 
there appears an infinite potential valley at x = 1. Transition in the direction — )■ 1 
becomes exponentially difficult with respective to ANfi. 
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Figure 4: Adaptive landscapes in bi-stable systems under the joint effects of drift, mu- 
tation, and two types of selection (shown respectively as solid lines in the two figures). 
In both cases ANfi = ANu = 0.1; the dotted line in each figure shows the landscape 
without selection as a reference, (a) The potential landscape of the first type (solid) is 



given in Eq.(30). Compared to the neutral potential (dotted), it skews in a way that 
shows the selectional preference of Ai. (b) The potential landscape of the second system 



is expressed in Eq.(32). The deeper valley (compared to the neutral landscape) shows the 
homozygotic advantage over heterozygotic types by selection. Bi-stability is more salient 
after introducing this advantage (s > 0). 
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approximation in Eq.(|34|) in the range 
p = 0.000003, 



Figure 5: Comparison between the numerical solution of Eq.(|33|) and the analytical 

< ANs 
100. When ANs < 1, 



ical solution accurately; the minor error for s 



< 1. Other parameter settings are: 
the leading term approximates the numer- 
= comes from the remaining terms of the 
expansion Eq.(24). It can certainly deal with s < 4/i. As ANs increases, the expansion 
of its exponential terms converges slower. So the error of the total expansion becomes 
major. 
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